/*
 * Copyright (c) 1997, 1998, 2003, 2006 Aleksandar Samardzic
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. The name of the author may not be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include "plane.h"

/* Plane_CalcCoeffs - calculate plane equation coefficients from three
 * points */
void
Plane_CalcCoeffs(this, p0, p1, p2)
     PPlane         *this;
     PVector         p0;
     PVector         p1;
     PVector         p2;
{
	PVector         normal,
	                P,
	                Q;
	double          mul,
	                den;

	/* calculate plane normal */
	Vector_Sub(P, p1, p0);
	Vector_Sub(Q, p2, p0);
	Vector_CrossProduct(normal, P, Q);
	Vector_Normalize(normal, mul, den);

	/* calculate coefficients */
	this->A = normal[X], this->B = normal[Y], this->C = normal[Z];
	this->d = Vector_DotProduct(normal, p0);

	/* determine largest projection axis */
	this->ax =
	    (fabs(this->A) >
	     fabs(this->B)) ? ((fabs(this->A) >
				fabs(this->C)) ? X : Z) : ((fabs(this->B) >
							    fabs(this->
								 C)) ? Y :
							   Z);
	this->ax++, this->ax %= 3, this->ay = (this->ax + 1) % 3;
}
